Intermingled basins in coupled Lorenz systems 



(N 

O 
(N 

O 

m 



in 



Sabrina Camargo ^, Ricardo L. Viana ^, and Celia Anteneodo -^'^ 

^ Department of Physics, PUC-Rio, Rio de Janeiro, Brazil 
^ Department of Physics, Federal University of Parana, Curitiba, Brazil 
^ National Institute of Science and Technology for Complex Systems, Rio de Janeiro, Brazil. 

(Dated: January 31, 2012) 

We consider a system of two identical linearly coupled Lorenz oscillators, presenting synchro- 
nization of chaotic motion for a specified range of the coupling strength. We verify the existence 
of global synchronization and antisynchronization attractors with intermingled basins of attraction, 
such that the basin of one attractor is riddled with holes belonging to the basin of the other attractor 
and vice versa. We investigated this phenomenon by verifying the fulfillment of the mathematical 
requirements for intermingled basins, and also obtained scaling laws that characterize quantitatively 
the riddling of both basins in this system. 
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I. INTRODUCTION 



The Lorenz system 

X = a{y — x), y = l3x — y — xz, z — —jz + xy, (1) 

for a — 10, (3 — 28, and 7 — 8/3, displays a chaotic 
attractor with the familiar butterfly-like shape [l| . It 
is often quoted as a paradigmatic system in nonlinear 
dynamics, since it displays many interesting dynamical 
properties of chaotic dissipative systems. Moreover its 
equations mimetize the dynamical behavior expected to 
occur in some physically relevant systems, as convection 
rolls in the atmosphere [3, single-mode lasers , and 
segmented disk dynamos Q. Coupled Lorenz systems 
could arise as well in the mathematical modeling of re- 
lated physical problems. The simplest case in the latter 
category is the coupling of two identical Lorenz systems. 

In identical coupled systems, even if chaotic, synchro- 
nization of trajectories may occur This phenomenon 
has been studied for more than two decades, motivating 
a wealth of analytical, numerical, and even experimen- 
tal results 5]. Synchronization of chaos, besides its own 
interest as a mathematical problem, finds applications 
for instance in secure communications The chaotic 
nature of the dynamics of one of the systems can be ex- 
ploited to code messages which could be sent to an iden- 
tical system through some form of coupling. If the latter 
system is synchronized with the former, the message can 
be securely uncoded. 

For two completely synchronized systems, cither peri- 
odic or chaotic, their dynamical variables are equal for all 
times. On the other hand, if instead of the difference, it 
is the sum of some of their dynamical variables that van- 
ishes, the two systems are said to antisynchronize. Due 
to phase-space symmetries, coupled Lorenz systems can 
exhibit both synchronized and antisynchronized states. 
Then, for secure communications purposes, the existence 
of another, antisynchronized, state is in principle a source 
of troubles since, depending on the initial condition, the 
receiver system could be tuned to the antisynchronized 



attractor. This situation can still be dramatically worsen 
when the riddling phenomenon occurs. 

As a matter of fact, multistable dynamical systems 
typically have a very complicated structure of basins 
of attraction, that may be delimited by fractal bound- 
aries Q- Suppose, for instance, that a dynamical system 
has two attractors, with the corresponding basins of at- 
traction sharing a common basin boundary in the phase 
space. If a ball centered at a given initial condition and 
with a radius equal to the uncertainty level intercepts the 
basin boundary, we cannot say a priori which attractor 
the system will asymptote to [1]. If that boundary is 
a curve, even if fractal, the final-state sensitivity prob- 
lem can be circumvented by decreasing the radius of the 
uncertainty ball (this can be done in experimental or nu- 
merical settings, by increasing the precision in determin- 
ing the initial condition in phase space). However, such 
reduction of uncertain initial conditions is not possible in 
the limit case in which the fractal boundary is area-filling, 
i.e., the (box-counting) dimension of the basin boundary 
gets close to the dimension of the phase space itself . 
In that limit case, the fraction of uncertain initial con- 
ditions will likely not decrease no matter how much we 
decrease the uncertainty balls of each initial condition. 
The latter situation occurs for riddled basins 0. 

From the mathematical point of view, riddled basins 
are observed in dynamical systems that exhibit an invari- 
ant smooth hypersurface with a chaotic attractor lying 
on it, another asymptotic final state, out of the invariant 
subspace, and negative Lyapunov exponent transverse to 
the invariant subspace with positive finite-time fluctua- 
tions [l0l - [l^ . Under the conditions above, riddfing orig- 
inates from the loss of transversal stability of unstable 
periodic orbits embedded in the chaotic attractor [l3| . 
despite the attractor being transversely stable in aver- 
age. In this context, attractors must be understood in 
the weak sense of Milnor [l^ . The transition associated 
to the flrst unstable orbit on the attractor that losses 
transversal stability determines the riddling bifurcation 
(see for instance Ref. [l5[ for an overview). Depending 
on the way these orbits loss stability and even on the 
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dynamics outside the invariant manifold, different bifur- 
cation scenarios and different forms of riddled basins can 
occur [13I [l6l - [l8l | (to cite a few examples). 

If riddled basins exist in a multistable chaotic system, 
their final states are utterly unpredictable, i.e. we can- 
not say - with any degree of certainty - which attractor 
the system will evolve to for long times [l^. The situa- 
tion, in this case, is similar to that for a random process, 
for which there can only be determined a probability for 
predicting the final state of the system. In fact, some 
phenomena formerly attributed to random variations in 
initial conditions can be also interpreted as a consequence 
of riddling [20|. 

The simplest case of riddling is when only one of the co- 
existing attractors have a riddled basin. However, when 
there is more than one invariant subspace, then more 
than one attractor can be riddled. In this case, the basin 
structure is called intermingled [9]]. 

The aim of the present work is precisely to show the 
existence of intermingled basins of attraction for the syn- 
chronized and antisynchronized states of two coupled 
Lorenz oscillators. In previous literature there are al- 
ready clues of such phenomenon. Kim and coworkers 
[2l| . in a work about anti-synchronization of coupled 
chaotic oscillators, point to the possibility of a riddled 
basin of synchronization in coupled Lorenz systems, but 
without going further on that issue. Furthermore, a one- 
dimensional reduction of the Lorenz system (to a piece- 
wise approximation to the well-known Lorenz map) was 
low-dimensional enough for an analytical treatment to 
be feasible and show the riddling of the synchroniza- 
tion basin "2^. The verification of the transversal sta- 
bility conditions through direct methods (i.e., by mak- 
ing a linear stability analysis of each invariant subspace) 
is quite difficult in two coupled Lorenz systems, since 
the phase space is six-dimensional. Then, we investigate 
those properties numerically. We also characterize quan- 
titatively the riddled basins by means of the scaling laws 
giving the probability of making wrong predictions on the 
final state of the system, with respect to two quantities 
of interest: (i) the phase-space distance to the invariant 
subspace; and (ii) the uncertainty radius for each initial 
condition [2^. We have verified that, for both quanti- 
ties, the probability scales as a power-law, as required 
for riddled basins. 

The rest of the paper is organized as follows: Sectionllll 
describes the coupled system of Lorenz oscillators, as well 
as the existence of both synchronized and antisynchro- 
nized states. Section IIIII presents a preliminary discus- 
sion of the basins of attraction of both the synchronized 
and antisynchronized states. The mathematical proper- 
ties required for riddled basins and the necessary tools are 
the object of Section |IVl Section |V] discusses the quanti- 
tative characterization of riddled basins through scaling 
laws and the theoretical results supporting them. The 
last Section contains our conclusions and final remarks. 



II. COUPLED LORENZ SYSTEMS 

Many different coupling schemes are possible for two 
identical Lorenz systems [23). We have chosen, for 
symmetry reasons, a diffusive coupling through the z- 
variable, as follows 



xi ^ a{yi - xi), 

yi = (3xi - 2/1 - xizi, 

zi= -"fzi + xiyi + s{z2 - Zl), 

X2 = a{y2 - X2), 

V2 = fiX2 - y2 - X2Z2, 

Z2= -JZ2+X2y2+£{zi- Z2), 



(2) 



where we will use the same values for a, /3, and 7, as in 
the uncoupled case, and e is the coupling strength. 

On considering the dynamical behavior of the coupled 
system, it is convenient to perform the changes of vari- 
ables 



X 



[X2 - Xi) 



{X2 + Xi) 



y = 

Y 



{yi - yi) 
2 

(j/2 +2/l) 



iz2 - Zl) 

2 

{Z2 + Zl) 



(3) 



2 ' 2 ' 2 

after which the coupled system ^ becomes 

X = a{y — x), 

y = I3x~y~ {Xz -f Zx), 

z = -{-1 + 2e)z + Xy + Yx, 

X = a{Y-X), 

Y = PX -Y - {XZ + xz), 

Z = -jZ + XY + xy. 



(4) 



Whenever more convenient to the analysis, we will refer 
either to the new or the old variables. 

From inspecting Eqs. Q there follows that the dynam- 
ics of the coupled system is invariant with respect to the 
transformation {x,y,z) — > {—x, —y, —z). Hence the con- 
ditions X ^ y ~ z — Q define an invariant subspace Ads ■ 
one initial condition that belongs to this subspace gener- 
ates a trajectory in phase space that remains in Ms for 
any time. This three-dimensional subspace defines the 
complete (or global) synchronization manifold character- 
ized by xi = X2, yi = y2, zi = Z2. 

The dynamics in the invariant subspace AAg^ described 
by the variables {X, Y, Z), is governed by the equations 
of the uncoupled Lorenz system, hence there is a chaotic 
attractor As (butterfiy-like shape) lying in A^g. 

Analogously, due to the symmetry {X, Y, z) 
{—X, —Y, —z), the states for which X ^ Y = z = define 
another invariant subspace A4a (anti-synchronization 
manifold), associated to the attractor Aa-, in which 
(x, y, Z) follows the dynamics of the uncoupled system, 
i.e. Aa is a Lorenz chaotic attractor in Ma- 
There are also other symmetries already present in the 
uncoupled Lorenz system. Notice in Eqs. ([2]) that either 
{xi,yi) {-xi,-yi) or (2:2, 1/2) ^ {-X2,-y2) lead, to 
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four-dimensional invariant subspaces, while both symme- 
tries together lead to a two-dimensional invariant sub- 
space with a saddle point at the origin. Finally, included 
in this two-dimensional subspace, the lines at zi = zi 
and z\ = —zi also represent invariant subsets. We did 
not find any other relevant attractor other than As and 
Aa, which are attractors of the dynamics in the respec- 
tive subspaces M.s and A4a, and can become attractors 
for the whole phase space depending on their transversal 
stability. 



III. BASINS OF ATTRACTION 

In dynamical systems with more than one attractor, 
the corresponding basins may have fractal boundaries 
and even more complicated structures like the Wada 
property [25| . Accordingly, in the coupled Lorenz system 
([2]), the two coexisting attractors representing synchro- 
nized and antisynchronized states are expected to have 
such complex basin boundary structure. 

Since the phase space of the coupled system is six- 
dimensional, the visualization of the basins of attraction 
depends on convenient phase space sections or projec- 
tions. Figure [T] shows a section of the basin of the anti- 
synchronization (synchronization) attractor Aa (-^s ), for 
different values of the coupling parameter. 

Each initial condition was integrated using a fourth- 
order Runge-Kutta scheme with fixed timestep 10"'^ and 
for a time t = 10^, after which we determined to what 
attractor the corresponding orbit has asymptoted t26|] . If 
an orbit has asymptoted to an antisynchronized (synchro- 
nized) state in Ada (A^s), its initial condition was painted 
black (white) . Hence the area painted black (white) is a 
numerical approximation of a section of the basin of at- 
tractor Aa {As)- We considered 10^ initial conditions 
with xi = X2 — yi — y2 — 1-0 while z\ and zi were 
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FIG. 1: Section at xi = 0:2 = J/i = j/2 = 1-0 of the basins of 
synchronization (white pixels) and antisynchronization (black 
pixels) attractors of the coupled Lorenz system, for e = (a) 
1.0, (b) 2.0, (c) 2.5 and (d) 2.8. 



randomly chosen in the interval [20, 24) according to a 
uniform probability distribution. 

For instance, for a coupling strength e = 1.0 [Fig. 
[IJa)], the section of the basin of attractor Aa is a se- 
ries of thin filaments stemming from the diagonal. The 
filaments are non-uniformly distributed and have a sug- 
gestive self-similar appearance. In fact, successive mag- 
nifications of [Fig. [IJa)] reveal similar patterns (see Ref. 
2l|). Such scenario is also observed for other values of e. 



as illustrated in Fig. Hfb-d) , even if some features change 
with e, such as the relative area of each basin, or the def- 
inition of the tongues anchored in the diagonal. Let us 
note that other cuts also display a tongue structure, as 
depicted in Fig. [5] for e = 2.0. 
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FIG. 2: Section at zi = Z2 = 22.0 and (a) 3/1=2/2 = 1.0, 
xi,X2 random in [—1,3), (b) xi = X2 = 1.0, y\,y2 random 
in [—1, 3) of the basins of synchronization (white pixels) and 
antisynchronization (black pixels) attractors of the coupled 
Lorenz system, for e = 2.0. 

The structure of the basins of attraction is indeed ex- 
pected to be altered by the coupling strength. As an 
example, in Fig. [3] we show that, for a given initial con- 
dition [xi = yi = zi = 1.0 and X2 = y2 = Z2 — 0.5) inte- 
grated up to time 10^, the trajectories in the subspace of 
each oscillator are distinct for e = 0.5, while for e = 1.0 
trajectories tend to coincide due to synchronization. In 
the latter case, the overlapping segments reproduce a cut 
of the familiar attractor of the single Lorenz system, since 
for synchronized orbits, the evolution proceeds towards 
the attractor in Ais which is defined hyx — y — z — 0, 
and X = xi — X2, Y = yi = y2, Z = zi = Z2 follow the 
dynamics of an uncoupled system, as described above. 
Differently, in the former case (e = 0.5), the trajecto- 
ries of each system depart from those of the uncoupled 
system. 

Moreover, the observation of synchronized or antisyn- 
chronized states depends on the coupling strength. Re- 
call that the existence of A4s (he. the synchronized state 
being a possible solution of the coupled equations) does 
not mean necessarily that synchronized states, and in 
particular states in its attractor As, can be observed 
in numerical simulations. This occurs only if there is 
transversal stability, in the sense that any infinitesimal 
displacement along directions transversal to M-s decays 
exponentially with time. 

Let us remark that, due to the symmetry of the 
equations with respect to synchronized/antisynchronized 
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FIG. 3: Trajectories of the coupled Lorenz system for the 
same initial condition {xi = ?/i = zi = 1.0 and X2 = y2 = 
Z2 = 0.5) up to t = 100 and different coupling values: (a) 
e = 0.5; (b) 1.0. In each case, the time evolution of the 
differences of coordinates are also shown (c)-(d) for e = 0.5 
and (e)-(f) for e = 1.0. 



states, comments for attractor As are also valid for Aa- 
In order to visualize the existence of a transversely sta- 
ble synchronization manifold, we consider the differences 
Xi — X2-, 2/1 — J/2, and zi — Z2, which must vanish if a 
synchronized attractor is achieved. For e > 0.7, zi — Z2 
vanishes [Fig. IHc)], while the other two differences may 
also vanish (global synchronization) or not (local syn- 
chronization) [Fig. lU^a) and (b)]. Similar plots are ob- 
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FIG. 



4: Difference between the coordinates (a) xi — X2; (b) 
y2', (c) zi — Z2 of two coupled Lorenz systems at time 



t = 10"* , as a function of the coupling strength. One hundred 
initial conditions were randomly chosen (as in Fig. [ij for each 
value of e (varied in steps of 0.01). 



tained for the sums Xi + X2-, yi + y2, indicating that the 
basins of the synchronized and antisynchronized states 
are complementary to each other. 
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FIG. 5: (a) Fraction of initial conditions (over a total of 10'') 
yielding trajectories asymptoting synchronized fs or antisyn- 
chronized fa states as a function of time, for e — 2.0. We 
also indicated the fraction of trajectories reaching either state 
(fs + fa) or none of them (1 — /s + /a), (b) Fraction of initial 
conditions not reaching these states for different values of the 
coupling strength. 

We did not find any relevant attractor for the coupled 
system other than As and Aa- Besides the symmetry 
considerations at the end of Sect. II, we performed the 
following numerical experiment: we considered the ini- 
tial conditions used to plot the sections in Fig. [1] and, 
for each time t we computed the fraction of initial con- 
ditions that go either to As or Aa [Fig. El^a)]. The sum 
of these fractions rapidly approaches 100% [Fig. [SKa)], 
meaning that the fraction of initial conditions that do 
not asymptote to them goes to zero (filled squares in Fig. 
[5fa)), suggesting the existence of only two attractors for 
the coupled system. This conclusion has been observed 
to hold for £ > 0.7 as illustrated in [Fig. EKb)]. Oth- 
erwise, neither synchronized nor antisynchronized states 
are approached, as illustrated in Fig.|31[c)-(d) for e = 0.5. 
(Therefore, basin diagrams as those shown in Fig. [T]will 
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be left blank). 



IV. RIDDLED AND INTERMINGLED BASINS 

The standard requirements for the existence of a rid- 
dled basin are the existence of (i) a smooth invariant sub- 
space (of lower dimension than the phase-space) contain- 
ing a chaotic attractor, (ii) another asymptotic final state 
(not necessarily chaotic) out of the invariant subspace, 
(iii) negativity of the Lyapunov exponents transverse to 
the invariant subspace with (iv) positive finite-time fluc- 
tuations [9l-[T]|. which are associated to the transversal 
stability properties of unstable periodic orbits (UPOs) 
embedded in the attractor. For two symmetrically inter- 
mingled basins, the requirements for mutual riddling can 
be summarized as follows: 

1. There are invariant manifolds Sg and Sa contained 
in the phase space H. 

2. The dynamics on each manifold Sg^a has a chaotic 
attractor. 

3. 5s, a are transversely stable, meaning that the 
largest transversal Lyapunov exponent A_l is nega- 
tive. 

4. Although weak stability holds in average (condition 
3), UPOs embedded in the chaotic attractor are 
transversely destabilized. 

In Section IIIII we showed that conditions 1 and 2 are 
fulfilled for the coupled Lorenz systems. There exists 
two (three-dimensional) manifolds, 5^,0 = Ms.a^ in the 
six-dimensional phase space. They are invariant since 
trajectories starting in each subspace will remain there 
forever. Because the dynamics in each subspace coincides 
with that of the uncoupled map, then, it will evolve to- 
wards the respective well known Lorenz attractor As,a 
lying in the corresponding manifold. 

Moreover, for each invariant subspace, there are out of 
three transversal directions. Condition 3 means that the 
transverse Lyapunov exponents of typical orbits lying in 
the invariant subspaces (Ma and A4s) are all negative. 
The point in parameter space where they become posi- 
tive defines the blowout bifurcation [llj. To investigate 
condition 3, it suffices to consider the largest transversal 
exponent, denoted as A_l = limt_5.oo Ax(xo,t) < 0, where 
xq is an initial condition on the basin of attraction of 
either Aa or As- 

We computed Lyapunov exponents using two differ- 
ent methods. The Lyapunov spectrum was obtained fol- 
lowing the algorithm of Wolf et al. ^7[, with a Gram- 
Schmidt normalization step of 0.1. We integrated Eqs. 
([2]) using initial conditions given by = yi = X2 = 
2/2 = 1-0 and zi, zi were randomly chosen in the interval 
[20, 24) from a uniform probability function, as in Fig. [T] 
These initial conditions lead to trajectories that asymp- 
tote to either As or Aa- As a matter of fact, this is not 



relevant since both attractors have the same Lyapunov 
spectrum. 

The second method we used, and which can be applied 
to obtain only the largest transversal exponent, is to con- 
sider the time evolution of an infinitesimal displacement 
along a direction transversal to the synchronized sub- 
space A^s, which is given by 0, 



= lim - In— 

t^oc t (5(0)' 



(5) 



where (5(t) = \J {SxY' + {^vY + {^^Y is the norm of the 
transverse displacement, whose evolution is given by the 
variational equations for {x,y, z), setting x — y — z — 0, 
i.e., 



Sx — ce{Sy — Sx), 

Sy — f3Sx — Sy — XSz — Z6x, 

6z = -{-f + 2e)Sz + XSy + YSx, 

X = a{Y -X), 

Y = PX-Y-XZ, 

Z = -jZ + XY. 



(6) 



Figure [S] shows (in gray symbols) the three largest 
(infinite-time) Lyapunov exponents as a function of the 
coupling strength e, obtained by means of the algorithm 
by Wolf et aL, and the largest transversal exponent given 
by (O is indicated by a thick black line. One of the ex- 
ponents is always zero, corresponding to displacements 
along the trajectory. The largest exponent is practically 
always equal to ^ 0.9 and corresponds to the chaotic dy- 
namics on As (Aa)- The third exponent is the largest 
transversal exponent which we focus our attention on, 
both methods being in good accord in the region of in- 
terest (as shown in Fig. 15]). For the chaotic attractors in 
both synchronization and antisynchronization manifolds. 




FIG. 6: The three largest Lyapunov exponents of the coupled 
Lorenz system as a function of the coupling strength. Gray 
symbols correspond to the algorithm by Wolf et al., whereas 
the thick black curve is the result of the variational equations 
©. Inset: the largest exponent for a wider range of £. 
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the largest transversal exponent vanishes, changing sign, 
at El « 0.714 ± 0.005 and £2 « 3.061 ± 0.005, defin- 
ing the critical points of the blowout bifurcation. (There 
is also another critical value for large e, as can be seen 
in the inset of Fig. |6l but we will restrict our analysis 
to the lower range only). The largest transversal expo- 
nent is negative for ei < e < £2, hence, in that interval, 
condition 3 for intermingled basins is fulfilled. However, 
while the invariant subspaces Ms,a are stable in aver- 
age, with negative transversal Lyapunov exponents, there 
may be particular unstable periodic orbits embedded in 
the chaotic attractors .4s, a that are also transversely un- 
stable, with positive largest transversal Lyapunov expo- 
nent [2^ (condition 4). When trajectories come close to 
these unstable orbits, they will be repelled from the vicin- 
ity of the attractor. This will be reflected in positive val- 
ues of the finite-time largest transversal Lyapunov expo- 
nent [Tol - [l^ . Then we numerically computed the finite- 
time largest transversal Lyapunov exponents A_L(xo,t), 
by means of Eq. ^ but at finite t. For large enough t, 
one recovers the infinite time exponent A^, which does 
not depend on Xq , for almost all initial conditions in the 
attractors As,ai in contrast to the finite-time ones that 
may depend on the initial condition. 

We quantify the contributions of the finite-time largest 
transversal exponent by obtaining a numerical approx- 
imation to the corresponding probability distribution 
function P(A_l(xo, t)). We considered a large number of 
points in Ms (with x = y = z = Q, X = Y — 1.0, and Z 
randomly chosen in the interval [20, 24)), and discarded a 
transient. These were the initial conditions used to com- 
pute the time-t largest transversal Lyapunov exponents. 
Alternatively, we generated a single long chaotic trajec- 
tory (after the transient has elapsed) and divided it into 
time-t segments, using then the ergodicity of the dynam- 
ics to ensure that the conditions are randomly chosen 
according to the natural measure of the attractor. The 
results were essentially the same. 




-1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 



i^(f=30) 

FIG. 7: Probability density functions of the time-30 largest 
transversal Lyapunov exponent for different values of e. Ini- 
tial conditions were taken as in Fig. [1] The full lines are 
Gaussian fits. 



Figure [7] shows probability distribution functions 
(when t — 30) for different values of the coupling 
strength. In all the considered cases the distribution is 
nearly Gaussian and presents positive tails. Then, we 
computed the positive fraction of finite-time exponents, 
ip{t) = P{X^{t))d\^{t) > 0. The positive fraction is 
plotted in Fig. |S] as a function of the coupling strength 
for different values of the time-i interval used to sample 
the finite-time exponents. 
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e 



FIG. 8: Positive fraction of the largest time-t transversal Lya- 
punov exponent as a function of the coupling strength, for 
different values of t. 

For e ^ £1, finite-time exponents soon become positive. 
This is in agreement with the positivity of the infinite- 
time exponent shown in Fig. [B] for this region, and also 
with the fact that trajectories do not approach the in- 
variant subspaces, but are soon repelled, as already noted 
when we tried to plot the basins in Section lllli which is 
not possible for that parameter range. The positive frac- 
tion drops rapidly to 50%, which corresponds to the case 
for which the infinite-time exponent vanishes (consistent 
with symmetric P{X±{t))), and then drops below 50%, 
when the infinite-time exponent is negative. For s ~ 1.4, 
the positive fraction is minimal. For larger values of £, it 
increases, crossing the 50% level again for £ ~ £2, where 
the infinite-time exponent is again zero at that point (see 
Fig. in]). After that, the positive fraction tends to 1 gently 
with £, yielding a positive infinite-time exponent. This 
smooth behavior, different from the abrupt one in the 
lower limit of the region of negative A^ , is consistent with 
the observation, for £ > £2 , of a basin structure reminis- 
cent of those in Fig. [U although the filaments from the 
diagonal are not neat. Even if trajectories are ultimately 
repelled, they can spend long time intervals close to each 
attractor. 

The fraction of positive finite-time exponents is non- 
null. However, for the range £1 < £ < £2, that fraction 
decreases with t, as expected because the distribution 
of finite-time exponents collapses towards a Dirac delta 
centered at Aj_ in the long time limit. The decay is expo- 
nential, the faster, the closer to the minimum at £ ~ 1.4. 
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Hence, the absence of an abrupt decay of the positive 
fraction indicates a nonvanishing fraction for finite times. 
Then, from this analysis, condition 4 cannot discarded for 
any range within the interval £i < e < 62- This suggests 
that at least one of UPOs should be transversely unstable 
in that interval. 
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FIG. 9: (Color online) (a) Largest transversal Lyapunov 
exponent, as a function of e, for the particular unstable peri- 
odic orbits (UPOs) embedded in the Lorenz chaotic attractor 
(up to period 5), indicated on the figure by means of the 
sequence of symbols A, B denoting the turns around each un- 
stable fixed points C"*" and C~ of the Lorenz system. The 
curve for typical chaotic trajectories in the attractor is also 
shown, (b) Magnification of panel (a), (c) Stabihty intervals 
for each UPO, in order of increasing stability at e = from 
top to bottom: stable (full segment), unstable (dotted). The 
vertical lines indicate ei and £2- The symbols delimiting the 
stability intervals correspond to ^ = 1 (full) and —1 (hollow). 



Then we inspected the transversal stability of those 
orbits along the lines of periodic orbits threshold theory 
[28j . Once localization in phase space and periods of low 
period UPOs are available in the literature for the Lorenz 
system [2^, we computed Floquet multipliers [1^ [30| . 
Namely, we integrated Eqs. (jH]), to obtain the matrix Q 



such that 6{t) = Q5(0), with r the time period of the 
orbit and 5{t) — {5x, Sy, Sz) the column vector of trans- 
verse deviations. The eigenvalue of Q, /i, with maximal 
modulus furnishes Aj_ = In |/.i|/T, for a particular periodic 
orbit. Fig. ini shows the behavior of as a function of e 
for particular UPOs, up to period 5. UPOs are labeled 
by means of the sequence of symbols A, B denoting the 
turns around each unstable fixed point C*^ and C~ of the 
Lorenz system. Symmetric orbits obtained by exchang- 
ing A-f-s-B or with cyclic symmetry were omitted. 

One observes that the lowest period orbit AB (period 
2) appears to be the first in destabilizing the vicinity of 
£2, hence defining a riddling bifurcation. Then, between 
this point and £2 riddling can occur. This interval cov- 
ers most of the range £1 < £ < £2, except for a very 
small interval in the vicinity of £1. However, note that 
orbits of the type A"B, with n = 1, 2, . . ., have a maximal 
transversal Lyapunov exponents that increases with n in 
the vicinity of £1, hence shrinking the remaining small 
region of stability around £ ~ 1 . To confirm whether this 
region of strong stability (with no transversely unstable 
orbits) actually disappears would require the analysis of 
higher period orbits, a hard task for this system, since the 
number of UPOs increases exponentially with the integer 
period. 

Near the blowout bifurcation at £1, the low-period 
UPOs (up to period 5) destabilize for coupling strength 
either weaker or stronger than the critical value, but 
close to it. Let us remark that, differently to the cou- 
pled Rossler system studied by Heagy et al. here 
the ordering of the exponents of the lowest period or- 
bits in the neighborhood of £1 is inverted with respect to 
the uncoupled case as depicted in Fig. [Q] This implies 
that paradoxically the most stable orbits in the attrac- 
tor are those responsible for the transversal destabiliza- 
tion in this parameter region. A similar inversion occurs 
on some domains of the parameter space of a system of 
symmetrically coupled Rossler oscillators (isi nl\ . This 
characteristic turns difficult the determination of the rid- 
dling bifurcation (first destabilized orbit) related to the 
blowout at £1, apparently triggered by higher period or- 
bits. 

Furthermore, our outcomes point to a different nature 
of the blowout bifurcations at £1 and £2. Nearby £1, 
UPOs destabilize in its vicinity. Moreover, for all the an- 
alyzed orbits, the multiplier /i crosses the circle |yLt| = 1 
along the real positive semi-axis (associated to a pitch- 
fork bifurcation). This is in contrast to the scenario at £2, 
where there are orbits destabilizing far from £2 and with 
multiplier /x either -1-1 or —1. In particular, the first orbit 
AB loses stability with /x = 1 . The differences are consis- 
tent with the picture given by finite-time exponents, for 
instance in connection with Fig. [51 where a much abrupt 
behavior of the positive fraction was encountered near 
£1. 

The intervals where riddling can occur are delimited on 
one side by the blowout bifurcation and on the other by 
the riddling bifurcation. In our case there are two of such 



intervals and they apparently overlap, such that at least 
one UPO has lost transversal stability in the full range 
El < e < £2, although this would have to be confirmed 
by the analysis of high period orbits, it is supported by 
the analysis of finite-time Lyapunov exponents. 



V. SCALING LAWS FOR RIDDLED BASINS 

In this Section, we will focus on the determination of 
the scaling properties of the basins, which provide a mea- 
sure of their structure [13] ■ Let us focus on the black fil- 
aments in Fig. [T{b), which belong to the basin of the 
antisynchronization attractor. They are anchored at the 
diagonal line, which is a cut of the synchronization man- 
ifold Ms, given hyx — y~z = and containing a 
chaotic attractor As, while the antisynchronization at- 
tractor Aa lies elsewhere. In Fig. ITOfa) we portrait a 
schematic picture of that structure. The filaments of the 
basin of Aa are tongues anchored at points of As, and 
the complement of the filament set belongs to the basin 
of As- If an initial condition starts within any of these 
narrow tongues, even if it is very close to As, the resulting 
trajectory will asymptote to the other attractor. 

The set of basin filaments for Aa is expected to be 
self-similar by quite general grounds. Once the riddling 
bifurcation occurs for a given periodic orbit, it also oc- 
curs for every preimage of this orbit, yielding a dense set 
of tongue-like sets anchored at the corresponding preim- 
ages on As ■ The tongue-like shape is a consequence of the 
nonlinear terms in the equations describing the transver- 
sal dynamics. The characteristic feature of riddling is 
that those tongues have widths that tend to zero as we 
approach As- Hence the basin of As always contains 
pieces of the basin of the other attractor, regardless the 
transversal distance to As , so forming a fine structure of 
basin filaments (the same applying to Aa)- 

This fine structure can be quantitatively characterized 
by the following numerical experiment ilh . i23i |: let us 
consider the invariant manifold atx — y = z = and, 
depart from that manifold, for instance, by increasing 
up to a distance I = 2\z\ = \zi — Z2\ [as depicted by the 
red line segment in Fig. [TOFa')]. Then we evaluate the 
fraction Vi of points in that segment that belongs to the 
basin of As- We obtained a numerical approximation of 
this fraction by considering a number of initial conditions 
X = y = 0, \z\ = 1/2, X = Y = 1.0, and Z randomly cho- 
sen in the interval [20,24). If the trajectories did not 
synchronize (within a small tolerance) up to a time such 
that transients have elapsed and stationarity holds (typ- 
ically, t — 10^), we consider that they asymptote to the 
antisynchronization attractor Aa and, accordingly, they 
do not belong to the basin of .4s . If the latter is riddled 
with tongues belonging to the basin of Aa, for any dis- 
tance I (no matter how small) there is always a nonzero 
value of Vi- This fraction tends to zero as / — )■ 0. The 
fraction of length belonging to the basin of Aa (fraction 
of trajectories that do not synchronize) can be written as 
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FIG. 10: (a) Schematic figure showing the structure of riddled 
basins near the invariant subspace that contains a chaotic at- 
tractor. (b) Fraction of trajectories Pi, that asymptote to the 
antisynchronized state as a function of the distance I — 2\z\ 
to the synchronized state, for diflferent values of the coupling 
strength e. The full lines are least squares fits, (c) for 
£ = 1.5 and diflferent orientations of the deviation /. 



Pi, ~ 1 — V; , and is expected to scale with / as a power law 
Pi,{l) ^ l^, where 77 > is a scaling exponent. We inte- 
grated several initial conditions at the same distance I to 
the synchronization subspace and computed the fraction 
Pi, of initial conditions that do not synchronize, repeating 
this procedure varying the distance /. The results shown 
in Fig. llOr b) confirm the existence of a power law for this 
fraction, for many values of the coupling strength. 

The results do not vary appreciably when one departs 
from the synchronization manifold in other directions 
other than z. In Fig. [TUf c') we plotted, for the same 
coupling strength (e — 1.5), the fraction for initial 
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conditions with 
and also for |a; 

dom values of ?/, 2 such that d 



z\ = 1/2, X = y = (open squares), 
1/2, y = z = (filled circles) and ran- 



y.2 _|_ y2t 



+ = 1/2 

(filled squares) , obtaining essentially the same scaling ex- 
ponent. This scaling behavior is observed within the in- 
terval (61,62), below £1, no trajectories synchronize as 
seen in the previous sections, above 62, one observes a 
synchronized fraction but it does not change with I. The 
dependence of the numerically determined scaling expo- 
nents 77 on the coupling strength is depicted in Fig. 1111 
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FIG. 11: Scaling exponent 77 for the fraction of trajectories 
that asymptote to the antisynchronized state obtained by a 
numerical experiment (filled circles). For comparison, the the- 
oretical values, given by 7] — \X±\/D are also plotted (open 
squares) . 

An analytical expression for the exponent 77 was de- 
rived by Ott and coworkers for a simple model (piecewise 
linear non-invertible map) [23| . Their theoretical predic- 
tion arises from a diffusion approximation for a biased 
random walk that mimics the fluctuations of finite-time 
largest transversal Lyapunov exponents A_l (t) . They ob- 
tain the law P^, ~ l^, with 77 = |A_l|/Z? where D is the 
diffusion coefficient. This diffusion approximation is ex- 
pected to be valid near the blowout bifurcation {X± ~ 0) 
of an attractor with a riddled basin. The authors conjec- 
ture that a similar diffusion approximation, hence a simi- 
lar relation involving parameters A_l and D, must rule the 
scaling relation in a large class of systems. The distribu- 
tions shown in Fig. [7] already display a Gaussian charac- 
ter, which improves with larger time-i interval, consistent 
with the probability distribution function of independent 
random innovations, that, by the central limit theorem, 
is Gaussian. Additionally, we plot in Fig. [T^ the variance 
^Ai(t) '^^ probability distribution functions for X±{t) 
as a function of time, for different values of the coupling 
strength. As a matter of fact, the variance decays with 
time towards zero following asymptotically a power-law 
with exponent —1, as required for a normal diffusion pro- 
cesses, so validating the stochastic approach of Ott et at. 
Accordingly, the diffusion coefficient D can be estimated 
from the numerical curves, following f^^j^-j ~ 2D/t. The 



estimates 77 = |A_l|/D are plotted in Fig. [Til together 
with the numerical values. Numerical and estimated val- 
ues are in very good agreement in the proximity of the 
critical values, as expected [sj]. For intermediate values 
(0.75 < e < 2) there is a discrepancy, and the numerical 
exponent remains close to one (linear behavior), as also 
observed for other systems with intermingled basins ^3 ■ 
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FIG. 12: Time decay of the variance of the finite-time largest 
transversal Lyapunov exponent for different values of the cou- 
pling strength e. The lines are least squares fits of the function 
f{t) = 2D/t to the numerical points, for large t, allowing to 
estimate the diffusion coefficient D. 

Another scaling law typical of riddled basins is related 
to the fraction of uncertain initial conditions, with re- 
spect to their final-state [Slj. We may regard riddled 
basins as an extreme case of fractal basins, for which 
there is final-state sensitivity and the uncertainty fraction 
scales as a power-law with the uncertainty level, whose 
exponent gives a measure of the extreme final-state sen- 
sitivity due to riddling. Consider again the points at 
X = 7/ = and \z\ — 1/2 drawn in the phase space por- 
trait in Fig. I13f a). as described earlier, and choose ran- 
domly an initial condition Xq on that region. Now choose 
randomly another initial condition Xq with uniform prob- 
ability within an interval of length 2^ and centered at Xq 
[Fig. [T^Ta)]. If both points belong to different basins, 
they can be referred to as ^-uncertain [33l.[33j. 

The fraction of ^-uncertain points, or uncertainty frac- 
tion, denoted by (jp), is the probability of making a mis- 
take when attempting to predict which basin the initial 
condition belongs to, given a measurement uncertainty 
^. This probability scales with the uncertainty level as a 
power law of the form (p) ^ where > depends on 
both Xq and I. Numerical results are shown in Fig. [T^lb). 
The stochastic model of Ott et al. predicts a power-law, 
with exponent given hy (jj — A^/(4I?A||), a prediction 
that agrees with our numerical results close to the crit- 
ical points. However, for intermediate values, while the 
stochastic model predicts small (though nonzero) values 
(0 < 0.28), our numerical results for (p, as illustrated in 
Fig- fTWbl , yield much smaller values (by a factor greater 
than ten). As a matter of fact, for riddled basins, the ex- 
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FIG. 13: (a) Schematic figure showing the numerical determi- 
nation of the uncertainty fraction, (b) Fraction of uncertain 
initial conditions as a function of the uncertainty level, for 
different values of the coupling strength. The solid lines are 
least squares fits. 



ponent (j) should be rigorously zero (i.e. there would be 
no way to decrease the uncertainty fraction by decreas- 
ing the uncertainty level). Indeed, our results [Fig. [T^ b)] 
support this scaling law, with numerically obtained ex- 
ponents close to zero. 



VI. CONCLUSIONS AND FINAL REMARKS 

Riddled basins for the synchronization attractor of cou- 
pled Lorenz systems have been previously suggested in 
the literature but without a detailed characterization. In 
this work we offer numerical evidence that, for a spec- 
ified range of the coupling parameter (ei < e < £2), 
coupled Lorenz systems exhibit symmetrically riddled 
basins of attraction for synchronized and antisynchro- 
nized states. Since there are only two symmetric attrac- 
tors, their basins are intermingled. We firstly showed 
that the mathematical conditions for the existence of rid- 
dled basins are fulfilled, with the help of properties of 
finite-time largest transversal Lyapunov exponents and 
of the largest transversal exponent for particular or- 



bits. This is important as furnishes the sources of local 
transversal instability of the attractor even if stable in 
average. In a second place, we verified the existence of 
two scaling laws characterizing quantitatively the degree 
of uncertainty related to the riddled basins. These nu- 
merical results were compared to an analytical prediction 
(the stochastic model [23|), yielding a good accord where 
expected. Beyond the characterization of the structure 
of a riddled basin, these scaling laws allow to quantify 
the limitations to improve the ability in determining the 
final state of the system by increasing the accuracy level. 

Let us remark that intermingling, in particular of sym- 
metric basins, has also been observed in other systems 
with either continuous (e.g., mechanical system |12| and 
coupled Rossler oscillators 11 71) or discrete time dynam- 
ics (coupled logistic maps [iMl)- In the latter case, the 
analysis of the lowest period orbit was enough to furnish 
the conditions for the occurrence of riddling in certain 
parameter region. In fact, as anticipated by the results 
presented in Fig. [9l deepening in that point of view may 
furnish precise information on the nature of the bifurca- 
tions triggering riddling, although this may be a difficult 
task for the present system. As other perspectives for fu- 
ture work on this system, let us also mention the plausible 
occurrence, beyond the blowout bifurcation, of two-state 
on-off intermittency (l^ for which there is some evidence 
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Finally, it can still be worthy to explore other re- 
gions of phase space, as well as other ranges (negative or 
large values) of the coupling parameter. 

In any case, for applications, multistability is already 
a source of troubles. Still worst, the existence of inter- 
mingled basins of attraction for the synchronized and an- 
tisynchronized chaotic states of this system jeopardizes 
the solution of the problem of ensuring a given final state, 
since the initial condition determination is always done 
within a certain uncertainty level. With riddled basins, 
any uncertainty level, however small, lead to complete 
indeterminacy of the future state of the system. Hence 
in this case we cannot use synchronization of chaos for 
any practical purpose, since we will always be haunted 
by the existence of the another, antisynchronized state, 
with a basin intermingled with the basin of the synchro- 
nized state. Of course, the same difficulties concern the 
predictability of natural phenomena modeled by coupled 
Lorenz systems. Therefore, the importance of detecting 
the regimes where riddling can occur in a dynamical sys- 
tem. 
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